Spectra Story

Author

Adam Kemberling

Published

October 9, 2024

Code
# Libraries
library(tidyverse)
library(gmRi)
library(rnaturalearth)
library(scales)
library(sf)
library(gt)
library(patchwork)
library(lmerTest)
library(emmeans)
library(merTools)
library(tidyquant)
library(ggeffects)
library(performance)
library(gtsummary)
library(gt)
library(sizeSpectra)
library(ggdist)
conflicted::conflict_prefer("select", "dplyr")
conflicted::conflict_prefer("filter", "dplyr")
conflicted::conflicts_prefer(lmerTest::lmer)

# Processing functions
source(here::here("Code/R/Processing_Functions.R"))

# ggplot theme
theme_set(
  theme_gmri(
    rect = element_rect(fill = "white", color = NA), 
    strip.text.y = element_text(angle  = 0),
    axis.text.x = element_text(size = 8),
    legend.position = "bottom"))



# vectors for factor levels
area_levels <- c("GoM", "GB", "SNE", "MAB")
area_levels_long <- c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")
# table to join for swapping shorthand for long-hand names
area_df <- data.frame(
  area = c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "All"),
  survey_area = c("SS", "GoM", "GB", "SNE", "MAB", "Northeast Shelf"),
  area_titles = c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"))


# Degree symbol
deg_c <- "\u00b0C"

Storycrafting - Northeast Shelf Spectra

I think we need to step back a bit and find the main story again–which is about how size structure of the NES fish community has changed over time, whether that varies by region, and what covariates are associated with those changes.
The Aug draft you developed really focused on comparing regression models to consider drivers based on the fish community size spectrum as viewed from the perspective of the Wigley species set vs. the full species set.
However, I didn’t see figures of those spectra in that draft (maybe there are somewhere else and I’ve lost track). The recent slide deck really focuses on seasonal differences in the size spectra and influence of temperature. - KMills

Starting from Square Uno

I’d like for us to start by rebuilding the story of how the size structure has changed over time and regional differences.
As such, I would like for us to look at plots of the size spectrum slope for (1) the all-species length spectrum and (2) the Wigley species biomass spectrum from annual, seasonal and regional perspectives.
I am thinking this would be two panel plots (all-spp length, Wigley-spp biomass) with 5 subpanels each (NES, GOM, GB, SNE, MAB) that are set up to show lines by seasonal and annual time steps, much like the plot on the right of slide 7 in your most recent slide deck…except with NES added and an annual line added. Keeping the significant trend lines in is helpful. This is the first step, so when you have that ready, please share it with the group. - KMills

Code
# Data used for Wigley estimates
wigley_in <- read_csv(here::here("Data/processed/wigley_species_trawl_data.csv"))
finfish_in <- read_csv(here::here("Data/processed/finfish_trawl_data.csv"))

Aside: Biomass Fraction in Each Community

If we need to make a decision about which community we want to use, here is the different biomass totals for each:

Code
# glimpse(wigley_in)

# Load raw
# General tidying only, removal of strata outside our study area, Spring and Fall only
# (removes inshore and rarely/inconsistently sampled strata etc.)
# shrimps removed
trawl_basic <- read_csv(here::here("Data/processed/trawl_clean_all_data.csv"))

# Biomass totals
raw_totals <- trawl_basic %>% 
  distinct(est_towdate, cruise6, station, stratum, tow, comname,  biomass_kg) %>%   pull(biomass_kg) %>% 
  sum()

wig_biomass_total <- wigley_in %>% 
  distinct(est_towdate, cruise6, station, stratum, tow, comname,  biomass_kg) %>%   pull(biomass_kg) %>% 
  sum()

ffish_biomass_total <- finfish_in %>% 
  distinct(est_towdate, cruise6, station, stratum, tow, comname,  biomass_kg) %>%   pull(biomass_kg) %>% 
  sum()


tibble(
  "Community" = c("All (including crustaceans)", "All finfish", "Wigley Species Subset"),
  "Biomass Total" = c(raw_totals, ffish_biomass_total, wig_biomass_total)) %>% 
  mutate(`% of Total` = (`Biomass Total` / raw_totals)*100) %>% 
  gt()
Community Biomass Total % of Total
All (including crustaceans) 3781548 100.00000
All finfish 3708046 98.05629
Wigley Species Subset 3539719 93.60501

Size Spectra Summary Tables

Code
# Summary tables
lspectra_mod_ffish %>% 
  gtsummary::tbl_regression()   %>% 
  add_global_p(keep = T) %>% 
  bold_p(t = 0.05)  %>%
  bold_labels() %>%
  italicize_levels() %>% 
  as_gt() %>% 
  tab_header(title = "Full Finfish Community - Length Spectra Trends")
Full Finfish Community - Length Spectra Trends

Characteristic

Beta

95% CI

1

p-value

est_year 0.00 0.00, 0.00 0.2
survey_area

0.005
    GoM
    GB 1.4 -5.5, 8.4 0.7
    SNE 6.8 -0.24, 14 0.058
    MAB -6.1 -13, 0.91 0.087
season

0.2
    Fall
    Spring 4.1 -2.9, 11 0.2
est_year * survey_area

0.005
    est_year * GB 0.00 0.00, 0.00 0.6
    est_year * SNE 0.00 -0.01, 0.00 0.040
    est_year * MAB 0.00 0.00, 0.01 0.12
est_year * season

0.2
    est_year * Spring 0.00 -0.01, 0.00 0.2
survey_area * season

0.7
    GB * Spring -1.2 -11, 8.7 0.8
    SNE * Spring -5.2 -15, 4.7 0.3
    MAB * Spring -4.5 -14, 5.5 0.4
est_year * survey_area * season

0.6
    est_year * GB * Spring 0.00 0.00, 0.01 0.8
    est_year * SNE * Spring 0.00 0.00, 0.01 0.3
    est_year * MAB * Spring 0.00 0.00, 0.01 0.3
1

CI = Confidence Interval

Code
# Summary tables
lspectra_mod_wig %>% 
  gtsummary::tbl_regression()   %>% 
  add_global_p(keep = T) %>% 
  bold_p(t = 0.05)  %>%
  bold_labels() %>%
  italicize_levels() %>% 
  as_gt() %>% 
  tab_header(title = "Wigley Community - Length Spectra Trends")
Wigley Community - Length Spectra Trends

Characteristic

Beta

95% CI

1

p-value

est_year 0.00 -0.01, 0.00 0.030
survey_area

0.3
    GoM
    GB -3.5 -10, 3.6 0.3
    SNE -6.4 -14, 0.69 0.077
    MAB -5.9 -13, 1.2 0.10
season

0.4
    Fall
    Spring 2.8 -4.2, 9.8 0.4
est_year * survey_area

0.3
    est_year * GB 0.00 0.00, 0.01 0.3
    est_year * SNE 0.00 0.00, 0.01 0.094
    est_year * MAB 0.00 0.00, 0.01 0.13
est_year * season

0.4
    est_year * Spring 0.00 0.00, 0.00 0.4
survey_area * season

0.002
    GB * Spring 3.7 -6.3, 14 0.5
    SNE * Spring 6.8 -3.2, 17 0.2
    MAB * Spring -12 -22, -1.6 0.023
est_year * survey_area * season

0.002
    est_year * GB * Spring 0.00 -0.01, 0.00 0.5
    est_year * SNE * Spring 0.00 -0.01, 0.00 0.2
    est_year * MAB * Spring 0.01 0.00, 0.01 0.018
1

CI = Confidence Interval

Code
# Summary tables
bmspectra_mod_wig %>% 
  gtsummary::tbl_regression()   %>% 
  add_global_p(keep = T) %>% 
  bold_p(t = 0.05)  %>%
  bold_labels() %>%
  italicize_levels() %>% 
  as_gt() %>% 
  tab_header(title = "Wigley Community - Mass Spectra Trends")
Wigley Community - Mass Spectra Trends

Characteristic

Beta

95% CI

1

p-value

est_year 0.00 0.00, 0.00 0.020
survey_area

0.3
    GoM
    GB -2.8 -6.2, 0.71 0.12
    SNE -2.6 -6.1, 0.93 0.15
    MAB -2.6 -6.1, 0.93 0.15
season

0.5
    Fall
    Spring 1.2 -2.3, 4.7 0.5
est_year * survey_area

0.4
    est_year * GB 0.00 0.00, 0.00 0.12
    est_year * SNE 0.00 0.00, 0.00 0.2
    est_year * MAB 0.00 0.00, 0.00 0.2
est_year * season

0.5
    est_year * Spring 0.00 0.00, 0.00 0.5
survey_area * season

<0.001
    GB * Spring 1.7 -3.2, 6.6 0.5
    SNE * Spring 1.3 -3.6, 6.2 0.6
    MAB * Spring -7.3 -12, -2.4 0.004
est_year * survey_area * season

<0.001
    est_year * GB * Spring 0.00 0.00, 0.00 0.5
    est_year * SNE * Spring 0.00 0.00, 0.00 0.7
    est_year * MAB * Spring 0.00 0.00, 0.01 0.003
1

CI = Confidence Interval

Part 2: Relationships to External Factors

The following plots and summaries will relate the Wigley species spectra to Temperature and Live Landings

Code
# # Model Dataset for Wigley Species Subset
wigley_bmspectra_model_df <- read_csv(here::here("Data/model_ready/wigley_community_bmspectra_mod.csv"))

# Use one for plots
wigley_bmspectra_model_df <- wigley_bmspectra_model_df %>% 
  mutate(survey_area = case_when(
    survey_area == "GoM" ~ "Gulf of Maine",
    survey_area == "GB" ~ "Georges Bank",
    survey_area == "SNE" ~ "Southern New England",
    survey_area == "MAB" ~ "Mid-Atlantic Bight"),
  survey_area = factor(survey_area, levels = area_levels_long),
  season = fct_rev(season))

Bottom Temperature Changes

Bottom temperatures have increased in all four regions. The distribution of seasonal bottom temperatures are such that Gulf of Maine experiences a much smaller difference in seasonal temperature swings than the other regions.

Code
temp1 <- wigley_bmspectra_model_df %>% 
  ggplot(aes(est_year, bot_temp, color = season)) +
  geom_point(size = 0.8, alpha = 0.5) +
  geom_ma(size = 1, n = 5, ma_fun = SMA, linetype = 1) +
  scale_color_gmri() + 
  facet_wrap(~survey_area, ncol = 1) +
  scale_y_continuous(labels = label_number(suffix = "\u00b0C")) +
  labs(y = "Bottom Temperature", x = "Year", color = "5-Year Smooth")


# Plot the distribution as a raincloud plot
temp2 <- wigley_bmspectra_model_df %>% 
  ggplot(aes(fct_rev(survey_area), bot_temp, color = season, fill = season)) +
   ggdist::stat_halfeye(
     alpha = 0.4,
     adjust = .5, 
     width = .6, 
    .width = 0, 
    justification = -.2, 
    point_colour = NA
  ) + 
  geom_boxplot(
    fill = "transparent",
    width = .15, 
    outlier.shape = NA
  ) +
  ## add dot plots from {ggdist} package
  ggdist::stat_dots(
    ## orientation to the left
    side = "left", size = 1,
    ## move geom to the left
    justification = 1.12, 
    ## adjust grouping (binning) of observations 
    binwidth = .25
  ) + 
  # geom_point(
  #   ## draw horizontal lines instead of points
  #   shape = 95,
  #   size = 10,
  #   alpha = .1) +
  scale_color_gmri() + 
  scale_fill_gmri() + 
  scale_y_continuous(labels = label_number(suffix = "\u00b0C")) +
  coord_flip() +
  theme(legend.position = "none") +
  labs(y = "Bottom Temperature", x = "", color = "Season")

temp1 | temp2

The overlap in bottom temperature distributions is very similar to the overlap in mass spectra exponents.

Code
wigley_b_dist_plot <- wigley_bmspectra_model_df %>% 
  ggplot(aes(fct_rev(survey_area), b, color = season, fill = season)) +
  ggdist::stat_halfeye(
     alpha = 0.4,
     adjust = .5, 
     width = .6, 
    .width = 0, 
    justification = -.2, 
    point_colour = NA) + 
  geom_boxplot(
    fill = "transparent",
    width = .15, 
    outlier.shape = NA) +
  ## add dot plots from {ggdist} package
  ggdist::stat_dots(
    ## orientation to the left
    side = "left", size = 1,
    ## move geom to the left
    justification = 1.12, 
    ## adjust grouping (binning) of observations 
    binwidth = .005) + 
  scale_color_gmri() + 
  scale_fill_gmri() + 
  coord_flip() +
  labs(
    y = "ISD Exponent", 
    x = "",
    title = "Individual Body Mass Distribution (1g - Max)")



# Plot next to bottom temperatures
wigley_b_dist_plot | temp2

Which makes me lean towards thinking that the size distribution is responsive to the ambient water temperatures, and perhaps less stationary and responsive to long-term trends in some area.

Code
wigley_bmspectra_model_df %>% 
  ggplot(aes(bot_temp, b)) +
  geom_point(aes(color = season, shape = survey_area)) +
  geom_smooth(method = "lm", color = "black") +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(title = "Broad Relationship to Ambient Bottom Temperature")

Total Landings

Code
wigley_bmspectra_model_df %>% 
  distinct(est_year, survey_area, total_weight_lb) %>% 
  ggplot(aes(est_year, total_weight_lb/1e6)) +
  geom_col(alpha = 0.5) +
  geom_ma(aes(color = "5-Year Smooth"), size = 1, n = 5, ma_fun = SMA, linetype = 1) +
  scale_color_manual(values = "orange") +
  facet_wrap(~survey_area, ncol = 1) +
  scale_y_continuous(labels = label_number(suffix = "M")) +
  labs(y = "Total Landings (live lb.)", x = "Year", color = "Moving Average")


Driver Model

The following model has been used to determine the relationship between landings and seasonal bottom temperature on b.

lmer(
  b ~ survey_area*season*scale(roll_temp) + survey_area * scale(log(total_weight_lb)) + (1|est_year), 
  data = wtb_model_df)
Code
#### Model 2: Temperature & Fishing ####


# Drop NA's
wtb_model_df <- drop_na(wigley_bmspectra_model_df, total_weight_lb, bot_temp, b) %>%
  rename(area = survey_area) %>% 
  left_join(area_df) %>% 
  # Do rolling BT within a season * region
  group_by(survey_area, season) %>%
  mutate(
    roll_temp = zoo::rollapply(bot_temp, 5, mean, na.rm = T, align = "right",  fill = NA)) %>% 
  ungroup() %>% 
  mutate(
    yr_num = as.numeric(est_year),
    yr_fac = factor(est_year),
    survey_area = factor(survey_area, levels = area_levels),
    season = factor(season, levels = c("Spring", "Fall")),
    landings = total_weight_lb,
    yr_seas = str_c(season, est_year))


# Make the model
mod2 <- lmer(
  b ~ survey_area*season*scale(roll_temp) + survey_area*scale(log(total_weight_lb)) + (1|est_year), 
  data = wtb_model_df)


# summary(mod2)
# check_model(mod2)

Drive Model Sumamry

Code
mod2 %>% 
  gtsummary::tbl_regression()   %>% 
  add_global_p(keep = T) %>% 
  bold_p(t = 0.05)  %>%
  bold_labels() %>%
  italicize_levels() %>% 
  as_gt() %>% 
  tab_header(title = "Wigley Community - Bodymass Distribution & Drivers")
Wigley Community - Bodymass Distribution & Drivers

Characteristic

Beta

95% CI

1

p-value

survey_area

0.030
    GoM
    GB -0.13 -0.29, 0.02 0.088
    SNE 0.07 -0.06, 0.19 0.3
    MAB 0.05 -0.05, 0.15 0.3
season

0.2
    Spring
    Fall 0.06 -0.04, 0.17 0.2
scale(roll_temp) -0.01 -0.13, 0.10 0.8
scale(log(total_weight_lb)) 0.08 0.00, 0.15 0.048
survey_area * season

<0.001
    GB * Fall 0.17 0.00, 0.35 0.051
    SNE * Fall -0.12 -0.29, 0.06 0.2
    MAB * Fall -0.22 -0.39, -0.06 0.009
survey_area * scale(roll_temp)

0.007
    GB * scale(roll_temp) -0.24 -0.40, -0.08 0.004
    SNE * scale(roll_temp) 0.02 -0.14, 0.18 0.8
    MAB * scale(roll_temp) -0.10 -0.27, 0.07 0.3
season * scale(roll_temp)

0.7
    Fall * scale(roll_temp) 0.03 -0.11, 0.16 0.7
survey_area * scale(log(total_weight_lb))

0.4
    GB * scale(log(total_weight_lb)) -0.06 -0.17, 0.06 0.3
    SNE * scale(log(total_weight_lb)) -0.05 -0.17, 0.07 0.4
    MAB * scale(log(total_weight_lb)) -0.07 -0.14, 0.01 0.094
survey_area * season * scale(roll_temp)

0.031
    GB * Fall * scale(roll_temp) 0.16 -0.03, 0.36 0.10
    SNE * Fall * scale(roll_temp) -0.14 -0.33, 0.06 0.2
    MAB * Fall * scale(roll_temp) 0.06 -0.13, 0.26 0.5
1

CI = Confidence Interval

Temperature Relationship

Code
# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    mod2, 
    terms = ~ roll_temp*survey_area*season))   %>% 
  mutate(
    survey_area = factor(group, levels = area_levels),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####
trend_df <- emtrends(
  object = mod2,
  ~survey_area * season,
  var = "roll_temp")


# Drop effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = mod2,
  specs =  ~ survey_area*season,
  var = "roll_temp") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))




# Limit temperature plotting range
actual_range <- wtb_model_df %>% 
  group_by(season, survey_area) %>% 
  summarise(min_temp = min(bot_temp)-2,
            max_temp = max(bot_temp)+2,
            .groups = "drop")



# Plot over observed data
mod2_preds %>% 
  left_join(actual_range) %>% 
  filter((x < min_temp) == F,
         (x > max_temp) == F) %>% 
  left_join(select(mod2_emtrend, survey_area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = survey_area), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wtb_model_df,
    aes(roll_temp, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(survey_area~., scales = "free") +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(y = "Exponent of ISD",
       title = "Wigley Species, Body Mass Spectra",
       x = "5-Year Rolling Average Bottom Temperature")

Landings Relationship

Code
# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    mod2, 
    terms = ~ total_weight_lb * survey_area * season))   %>% 
  mutate(
    survey_area = factor(group, levels = area_levels),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####
trend_df <- emtrends(
  object = mod2,
  ~survey_area * season,
  var = "total_weight_lb")
#summary(trend_df, infer= c(T,T))


# Drop effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = mod2,
  specs =  ~ survey_area*season,
  var = "log(total_weight_lb)") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T),
    group = str_c(survey_area, "X", season))




# Plot over observed data
mod2_preds %>% 
  left_join(select(mod2_emtrend, survey_area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wtb_model_df,
    aes(total_weight_lb, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(survey_area~., scales = "free") +
  scale_color_gmri() +
  scale_x_continuous(transform = "log10", labels = label_log(base = 10)) +
  labs(y = "Body Mass Spectra Slope (b)",
       title = "Wigley Species, Body Mass Spectra",
       x = "Total Landings (lb.)")

Side Stories

Everything below is digging into different nuances of the approach or the community composition. These things are of secondary importance and will be moved out.

Spiny Dogfish Sensitivity

Spiny dogfish constitute a large majority fraction of biomass in the MAB during the Spring survey. Their numbers have also been increasing over time. This acts to shallow the spectra as they are one of the larger species routinely sampled by the trawl survey.

If we remove them from the estimation and re-run the slopes this is what changes.

How much can one species shift the spectra, and is this a definitive/conclusive proof that its the dogfish shallowing the slope in MAB?

Code
# # Run MAB with and without spiny dogfish
# # Run the year, season, area fits for the filtered data
# no_dogs_bmspectra <- group_binspecies_spectra(
#   ss_input = filter(wigley_in, comname != "spiny dogfish"),
#   grouping_vars = c("est_year", "season", "survey_area"),
#   abundance_vals = "numlen_adj",
#   length_vals = "length_cm",
#   use_weight = TRUE,
#   isd_xmin = 1,
#   global_min = TRUE,
#   isd_xmax = NULL,
#   global_max = FALSE,
#   bin_width = 1,
#   vdiff = 2)

# # Save those
# write_csv(
#   no_dogs_bmspectra,
#   here::here("Data/processed/wigley_species_nodogfish_bmspectra.csv"))
Code
# Read them in, do some plotting
no_dogs_bmspectra <- read_csv(
  here::here("Data/processed/wigley_species_nodogfish_bmspectra.csv")) %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels))


# Format a little
no_dogs_bmspectra <- no_dogs_bmspectra %>% 
  mutate(est_year = as.numeric(est_year))


# Join them together
dfish_results <- bind_rows(
   list(
     "Wigley Full" = wigley_bmspectra,
     "Wigley w/o Spiny Dogfish" =  no_dogs_bmspectra), 
   .id = "community") %>% 
  mutate(metric = "bodymass_spectra") %>% 
  mutate(
    survey_area = fct_relevel(survey_area, area_levels))




# Get trends for no dogfish
nodfish_mod <- lmer(
  formula = b ~ est_year * survey_area * season + (1|est_year),
  data = no_dogs_bmspectra)



dogfish_sens_predictions <- bind_rows(list(
  "bodymass_spectra-Wigley Full" = get_preds_and_trendsignif(bmspectra_mod_wig),
  "bodymass_spectra-Wigley w/o Spiny Dogfish" = get_preds_and_trendsignif(nodfish_mod)), 
  .id = "model_id") %>% 
  separate(model_id, into = c("metric", "community"), sep = "-") %>% 
  mutate(metric = fct_rev(metric))



# Plot the differences
ggplot() +
  geom_ribbon(
    data = filter(dogfish_sens_predictions, non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    data = dfish_results,
    aes(est_year, b, color = season),
    size = 0.8, alpha = 0.5) +
  geom_line(
    data = filter(dogfish_sens_predictions, non_zero == TRUE),
    aes(x, predicted, color = season),
    linewidth = 1, alpha = 1) +
  scale_fill_gmri() +
  scale_color_gmri() +
  facet_grid(survey_area ~ metric*community, scales = "free") +
  labs(
    x = "Year",
    y = "Exponent of Size Spectra",
    title = "Shift in Spectra Trend(s) with Exclusion of Spiny Dogfish")

Story Thoughts

Below is a table of some papers operating in the same area or around the same topic that are particularly relevant:

Code
shelf_papers <- tribble(
  ~"Author", ~"Year", ~"Conjecture",
  "Friedland et al.", 2024, "Body size has declined despite more abundance & biomass. Competition is preventing larger sizes.",
  "Krumsick", 2020, "Empirical slopes from trawl survey in labrador sea were steeper than theoretical from nitrogen stable isotope ratios"
)


shelf_papers %>% gt()
Author Year Conjecture
Friedland et al. 2024 Body size has declined despite more abundance & biomass. Competition is preventing larger sizes.
Krumsick 2020 Empirical slopes from trawl survey in labrador sea were steeper than theoretical from nitrogen stable isotope ratios

There are a couple core ideas about what changes have occurred or have been ocurring over the Northeast Shelf with relevance to the community that is sampled by the trawl survey: